%4-sep-2007. This code reads Julia and ACe data sets
% pick days with cont. common measurements
% 

%load c:\manoj\projects\plasma\HUA_PIU1.mat
load c:\manoj\projects\plasma\Julia_PLP.mat Julia_W fday w;
load c:\manoj\projects\ace\OMNI_ELEC ace_all ace_gsm_bz; % JUST ACE IEF
%load c:\manoj\projects\ace\OMNI_ELEC ace_gsm_bz; % JUST ACE IEF
load c:\manoj\projects\ace\pcn; mean_delta = 6.9444e-004;
%load c:\manoj\projects\ace\dmsp_pcp;PN_fday = f13_time; PN = f13_pcp;mean_delta = 0.0353652263;
pcn_fday = floor(PN_fday(:,1));
load c:\manoj\geomag\indices\aplist.mat;
%%

%Season Fliter
L = season_filter(PN_fday,11,1,3,1);
pcn_fday(~L) = [];
PN_fday(~L) = [];
PN(~L) = [];
%season filter over---------


% Hd = chebi11_filter();
N_seg = 1;
N_data = 0;
mjd_date = datenum(2000,1,1);

for i = 1: 749,
    
    if sum(diff(Julia_W(i).k))+1== length(Julia_W(i).k) &...
        sum(isnan(w(i,:))) <= 60, 
    L = pcn_fday == Julia_W(i).fday;
    if sum(L) > 0,
        
    pcn_time = PN_fday(L,1);%+(4*5)/(60*24);%Try to advance 20 PCN time minutes
    pcn_data = PN(L);
    
%
    L = isnan(pcn_data);
    if sum(L) < 10 & sum(L) ~= 0,
        pcn_data = interp1(pcn_time(~L),pcn_data(~L),pcn_time);
    end;
%

    fprintf('Number of missiing points on fday %d = %d\n', Julia_W(i).fday, sum(L));
   %if abs(mean(diff(ace_time))-0.0028) <= 2.7778e-005, % Use this with Burke data
   %if abs(mean(diff(ace_time))-0.0035) <= 2.7778e-005, %Use this with ACE min averages
   if abs(mean(diff(pcn_time)) - mean_delta) <= 2.7778e-005,% Use this with Thule 1 min PCN data
        JULI = w(i,Julia_W(i).k);
        PCND = interp1(pcn_time, pcn_data, fday(i,Julia_W(i).k));
%         JULI1 = filter(Hd,JULI);
%         PCND1 = filter(Hd,PCND);
       
%         subplot(311);
        L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
        mean_ap = mean(ap(L));
        
        L = ace_all(:,1) >= fday(i,Julia_W(i).k(1)) & ace_all(:,1) <= fday(i,Julia_W(i).k(end)); %IMF Bz data coems from ACE
        sel_data = ace_gsm_bz(L); % FOLLOWING THREE LINES ARE USED TO COMPUTE MEAN (withous using nanmean)
        L = isnan(sel_data);
        mean_imf_bz = mean(sel_data(~L));
        

        if sum(isnan(PCND)) < 1 & mean_ap >= 20,
            JULI_SEG(N_seg:N_seg+71) = JULI(1:72);
            PCN_SEG(N_seg:N_seg+71) = PCND(1:72);
            TIME_SEG(N_seg:N_seg+71) = fday(i,Julia_W(i).k(1:72));
            N_seg = N_seg+72;
            N_data = N_data+1;
            data_length(N_data) = length(JULI);
            AP_SEG(N_data) = mean_ap;
            IMF_SEG(N_data) = mean_imf_bz;
               
         end;
   else,
       fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
    end;
    end;
end;
end;

%[Cxx,F] = mscohere(JULI_SEG,ACE_SEG,hanning(72),0,72,0.0033333);
[Cxx,F] = mscohere(JULI_SEG,PCN_SEG,hanning(72),0,72,0.0033333);
 
figure1 = figure;

axes('Parent',figure1,'XTick',[10 20 30 50 70 100 200 300 500 1000],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);

box('on');
hold('all');
xlabel('period in minutes');
ylabel('coherence');
hold on;
semilogx((1./(60*F)),Cxx,'r','LineWidth',2);
grid on;
alpha = 0.95;
ci = 1- (1-alpha)^(1/(N_data-1));
aa = axis;
h=line([aa(1) aa(2)],[ci,ci],'LineStyle','-.','color','r');
axis([aa(1) aa(2) 0 1]);
